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INTRODUCTION 


In  this  paper,  we  consider  the  following  type  of  harvesting 


problem.  An  animal  population  is  divided  into  two  stocks:  an 


underlying"  population  and  a "surface"  population.  We  assume 


that  there  is  a natural  exchange  between  the  two  population  levels 
The  predator  or  harvestor  affects  only  the  "surface"  population 
and  does  not  influence  the  "underlying"  population  directly. 


Such  a situation  occurs 


Tropical  Tuna  Fishery  (e.g.,  IATTC , 1975).  In  this  case,  tuna 


associate  with  porpoise  schoools 


those  tuna  associated  with  porpoise.  Consequently,  the  underlying 
population  of  tuna  is  not  sampled  by  the  fishery.  One  may  wonder 


what  information  measurements  on  the  surface,  harvested  population 


provides  about  the  unobservable  underlying  population.  Furthermore 


it  is  interesting  and  important  to  know  if  the  standard,  linear 


population 


Clark  and  Mangel  (1977)  studied  some  of  these  questions 


as  they  relate  to  the  tuna  fishery.  They  constructed  a number  of 


of  the  models  exhibited  a bifurcation 


Namely,  the  steady  state  tuna  level,  N 


characterises  natural  association  and  growth  rates 


while  the  second,  E , characterizes  fishing  effort.  In  some  of  the 


models,  no  bifurcation  occurred.  In  other  models,  if  a is  less 


than  a critical  value  ac  , Clark  and  Mangel  showed  that  there  is 
a non-trivial  steady  state  N , E)  for  all  E.  When  a > ac, 
there  exists  a non-trivial  steady  state  only  if  E < E^,  some 
bifurcation  value  of  effort.  If  E > Eb  , then  = 0 is  the  only 

steady  state.  Clark  and  Mangel  used  a steady  state  analysis.  Here 
we  will  study  the  dynamical  equations  in  more  detail,  and  thus 
extend  the  previous  analysis. 

When  the  only  steady  state  is  the  trivial  one,  the  population 
will  approach  extinction  as  time  increases.  However,  for  many 
situations,  the  time  to  reach  N(t)  =0  is  not  as  relevant  as  the 
time  to  reach  N(t)  =5  , where  6 is  a given  (low)  population 
level.  This  is  true  for  the  following  reasons.  First,  the  predator 
will  cease  to  seek  the  prey  once  the  prey  becomes  scarce  enough. 

Second,  as  N(t)  -*  0 a model  in  which  N is  a continuous  variable 
is  no  longer  valid,  since  the  small  value  of  N demands  a 

discrete  model.  Third,  if  the  population  level  becomes  very  low, 
the  entire  dynamics  of  the  species  may  shift  into  some  other, 
stable  configuration  (e.g.,  a second  predator  may  become  important, 

(see  Holling,  1973) . We  introduce  a time  T^ (Nq,  E)  : 

i 

T.(N  , E ) = min[t:N(o)  = N ,N(t)  £ 6,N(S)>6,for  S<  t.E(t)=E]  (1) 

6 o o 


Then  T (Nq,E)  is  the  first  time  that  the  population  goes  below 

the  level  6 , given  that  it  starts  at  level  N and  is  harvested 
3 o 


with  constant  effort  intensity  E.  In  this  paper,  we  will  provide 


techniques  for  the  approximate  calculation  of  T^(Nq,  E)  . 

All  predator-prey  systems,  of  course,  operate  in  a random 
environment.  Clark  and  Mangel  did  not  consider  any  stochastic 
effects.  When  stochastic  effects  are  included,  the  proper  des- 
cription of  the  system  involves  probabilities  of  certain  events. 

Instead  of  T<$(No,  E)  we  must  consider  the  ensemble  average 

T6(No,  E)  = <T6(Nq,  E) > (2) 

We  will  show  how  T£(No,  E)  can  be  calculated. 

In  most  biological  systems,  the  modeling  of  stochastic  effects 
is  a difficult  problem.  We  present  a phenomenological  approach, 
in  which  deterministic  equations  are  reinterpreted  in  a birth  and 
death  framework.  This  technique  provides  a way  to  model  stochastic 
effects.  A second  approach  is  outlined  in  the  appendix. 

In  § 2 , we  rederive  model  B of  Clark  and  Mangel  (1977)  in 
a more  general  setting.  We  non-dimensionalize  the  equations,  so 
that  the  appropriate  parameters  become  clear.  The  validity  of 
the  steady  state  assumption  of  Clark  and  Mangel  can  then  be  studied. 

We  also  discuss  the  (N,  Q)  phase  plane,  where  N is  the  under- 
lying population  level  and  Q is  the  surface  population  level. 

Next,  we  show  how  to  incorporate  stochastic  effects  into  the  model 
and  formulate  the  probabilistic  questions  of  interest. 


3 


In  §3  , we  construct  solutions  of  the  deterministic  kinetic 
equations  for  N,Q  by  using  a singular  perturbation  technique 
(Lin  and  Segel,  1973).  In  particular,  we  give  a simple  technique 
for  the  calculation  of  Tg(NQ,  E) , accurate  to  order  1/E  . In 
§4  , we  consider  stochastic  effects  in  the  steady  state  approxi- 
mation Q = 0 . Exact  solutions  of  many  problems  are  possible. 

In  particular,  we  can  calculate  (Nq , E)  exactly.  Finally, 

in  §5  , we  discuss  stochastic  problems  in  the  (N,  Q)  phase 
plane.  The  relevant  theory  for  the  stochastic  problem  is  given  by 
Ludwig  (1975)  and  Mangel  (1977).  We  will  sketch  the  applica- 
tion of  those  theories,  but  will  not  go  into  great  detail.  There 
are  two  appendices.  In  the  first,  we  prove  a theorem  that  generalizes 
the  work  of  Clark  and  Mangel.  In  the  second,  we  provide  an  alterna- 
tive technique  for  the  modeling  of  fluctuations. 
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SECTION  2 

DERIVATIONS  AND  SCALING,  STOCHASTIC  MODELS 

We  generalize  the  model  of  Clark  and  Mangel,  without  the 
steady  state  assumption;  and  then  non-dimensionalize . This 
procedure  will  allow  us  to  study  the  validity  of  the  steady  state 
assumption.  Next,  we  introduce  the  stochastic  model  and  certain 
associated  probabilistic  questions. 

THE  AGGREGATION  MODEL  AND  SCALING 

Let  N (t ) , Q(t)  denote  the  underlying  population  level  and 
surface  population  level,  respectively,  at  time  t-  The  funda- 
mental assumption  is  that  underlying  population  aggregates  at 
the  surface  at  a rate 

Ra  = ctN  (1  - Q/Q)  ( 

In  equation  (3) , a is  a parameter  with  units  (time)  ^ J it  is  the 
association  rate.  Q is  a constant.  One  might  view  Q as  the 
carrying  capacity  of  the  surface  level. The  important  point  about 
the  form  of  (3)  is  that  the  surface  population  level  is  limited, 
regardless  of  the  value  of  N (Clark  and  Mangel,  1977).  In  fact,  the 
qualitative  results  reported  here  will  hold  for  any  aggregation 
model  in  which  the  surface  population  level  is  limited  (see 
appendix  A) . 

We  assume  that  the  population  has  a natural  growth  rate 
gN(N)  . In  many  instances,  we  assume  that  gN  is  logistic: 


. 


gN(N)  = rN (1  - N/N)  . 


We  also  assume  that  the  surface  population  is  harvested  at  a rate 
bE  , whera  b is  a constant  and  E is  fishing  effort.  We  let 
Xq  denote  the  average  fraction  of  the  surface  level  removed 
during  one  harvesting  operation.  The  dynamical  equations  for  sur- 
face and  underlying  stock  evolution  are: 

dN  _ „ ~M/,  Q . 
dt  gN  " aN^1  " - ^ 


= otN  ( 1 - - ) - bEX  Q . 
dt  Q ° 


We  now  show  that  these  lead  to  the  results  of  Clark  and  Mangel. 
We  assume  that  dQ/dt  = 0,  i.e.,  the  surface  level  rapidly 
achieves  a steady  state.  Then 


(4) 


(5) 


(6) 


ctNQ 


ss 


QbEX  + «N 
o 


(7) 


The  rate  at  which  the  surface  population  is  harvested  is: 


bEX  aNQ 

Y = bEX  O • K = - 

° ctN  + bEX_CT 


Equation  (8)  is  the  result  (15B)  of  Clark  and  Mangel.  The  steady 
state  assumption  assumes  that  deviations  from  Q die  off  on  a 

S 5 


(8) 


I 


rapid  time  scale.  In  order  to  calculate  this  time  scale,  we  non- 
dimensionalize  the  equations. 

Let 


x = rt 


Then  equations  (5,  6)  become 


H = n(l  - n)  - an(l  - q) 


(J\*Z 

\agjdx 


= n (1  - q)  - tt  q 


l 


We  define  aS  = 1/e  and  T/8  = X so  that  (11)  becomes 


e = n ( 1 - q)  - Aq. 


When  ag  = — — is  large,  e is  small.  For  the  tuna  fishery,  we 
rQ 

believe  that  e is  small  (Clark  and  Mangel,  1977). 

When  e = 0,  the  steady  state  q satisfies 


q = 


X + n 


(12a) 


Let  ^ denote  the  deviation  of  q (x  ) from  q : q (x  ) = 


_ A 

q + q (t  ) 


i 


fl 


Then, 


A 

e ^ = 

E dx 


-q (n  + A) 


Equation  (13)  allows  the  calculation  of  the  rate  of  relaxation  of 
deviations  from  q , (see  below) . 

THE  (n,  q)  PHASE  PLANE  AND  BIFURCATION  PICTURE 

We  now  consider  the  phase  plane  described  by  equations  (10) 

and  (12).  Steady  states  (or  equilibria)  are  determined  by  setting 

the  left-hand  sides  equal  to  zero.  We  obtain: 

n 

q A+n 

0 = n(l  - n)  - 


The  origin  (0,  0)  is  always  a steady  state.  The  other  steady 
states  are  the  roots  of : 


0 = (1  - n) (A  + n)  - aA 


0 = -n  + n(l  - A)  - aA  + A 


Hence,  there  will  be  another  real  steady  state  if 


D = (1  - A)“  - 4 A (a  - 1)  > 0 


We  distinguish  two  cases.  In  the  first  case,  a < 1 . Then 


D > 0 for  all  values  of  A . Hence,  we  always  observe  another 
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steady  state.  It  is  easy  to  show  that  there  is  only  one  non- 
negative steady  state.  The  phase  plane  is  sketched  in  figure  1. 
When  a > 1 , a more  interesting  phenomenon  occurs.  We  obtain  two 
non-trivial  steady  states  if 


a < 1 + 


(X  - 1) 
4X 


(18) 


The  phase  plane  in  this  case  is  shown  in  figure  2.  The  condition 
D = 0 allows  us  to  construct  a "bifurcation"  diagram  in  the  (a,  X) 
plane.  Such  a diagram  is  shown  in  figure  3.  We  note  that  this 
bifurcation  diagram  corresponds  to  the  "fold"  catastrophe  and  not 
the  "cusp"  catastrophe  (compare,  e.g.,  Jones  and  Walters  (1976)). 

The  bifurcation  value  of  X is  determined  by  D = 0,  i.e., 


a=l  + 


(xB  - 1) 


4X 


(18a) 


B 


STOCHASTIC  MODELS 

At  the  present  time,  there  is  much  discussion  concerning 
the  "proper"  way  of  introducing  stochastic  models  into  biological 
problems  (e.g.,  Turelli  (1977)  or  Ludwig  (1975)).  Ideally,  one 
would  start  with  a discrete  model  and  take  appropriate  limits  to 
obtain  deterministic  and/or  stochastic  equations  (Ludwig,  (1975) . 

Here,  we  shall  take  a different  approach  and  will  formally  re- 
interpret the  kinetic  equations  (10,  12)  to  obtain  a stochastic  model, 
as  is  done  by  Mangel  (1977) . 


I 


In  (19,  20),  we  have  rescaled  N,  Q so  that  they  represent  population 
level  in  numbers  (rather  than  biomass).  Equations  (19,  20)  are 
assumed  to  represent  the  evolution  of  the  average  value  of  random 
variables  N(t),  Q(t)  that  satisfy  stochastic  kinetic  equations 

(21) 


dN  = bN  (N , Q)  dt  + 4^  dWj 
d£f  = bQ(N,  Q)  dt  + ^ dW. 


In  (21,  22) , we  use  the  convention  of  summation  over  repeated 
indices.  The  process  W(t)  is  a Wiener  process.  Although  we  are 
using  an  Ito  equation,  we  shall  use  the  Stratonovich  inter- 
pretation (Wong  (1971)  , (Turelli  (1977)).  The  functions  A1-1  (N,  Q) 


are  defined  bv. 


i 


... 


i 1 


L I 


where  x\  X^  , Xk  are  ft  or  ft  . In  order  to  calculate  these 
variances,  we  reinterpret  (19,  20)  as  follows: 


,N , N ,N 
b (N,  Q)  = b+  - b_ 


(24) 


b°(N,  Q)  = bj  - b° 


(25) 


where 


bj  = rN  + ^ 


,N  _ rN' 


+ aN 


(26) 


N 


. Q _ ~KT 

b+  - aN 


bQ  _ aNQ  [ 
Q 


bEXQQ 


(27) 


We  view  b^  , b+  as"birth  rates";  b^  , b^  as"death  rates". 


/v  /v 


Consider  the  increments  AN  = N(t+At)  - N(t),  AQ  = Q(t+At)  - Q(t), 
given  that  N (t)  = N , Q(t)  = Q , we  assume  that 


Pr  |a&  = 1,  AQ  = o|  = rNAt  + o(At) 


Pr 


|Aft  = 1,  AQ  = -l|  = At  + o ( At) 


2 

Pr  { Aft  = -1,  AQ  = 0 } = At  + o (At) 

, 1 ’ N 

Pr  | Aft  = -1,  AQ  = l}  = aNAt  + o(At) 

Pr  | Aft  = 0,  AQ  = -l|  = bEXQQAt  + o(At) 


(27a) 


{ 


Pr  {all  other  transitions j = o(At) 
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We  now  consider  a number  of  probabilistic  quantities  connected 
with  the  stochastic  equations.  We  use  the  notation,  b1  for  bN 

or  b^  and  a^  for  the  elements  of  a . We  denote  by  x^  = N or  n 

2 12 
and  x =$  or  q. First  consider  the  density  for  (x  , x ) : 


v(x1,x2,t)dx'1'dx2  = Pr {x1  (t)£  (x1  ,x1+dx1)  , x"  (t)€(x2 ,X2+dx2 ) } (34) 


Then  v satisfies  the  forward  equation  (Ludwig,  1975) : 


vt  = ' ((bi  + ci)v)i 


In  equation  (35) , 

i / \ Id  i j 

c (x)  = j — T a 

* 3x: 

The  function  c(x)  arises  in  the  Stratonovich  interpretation  of 
the  equations  (31,  32).  In  (35),  subscripts  indicate  partial  dif- 
ferentiation and  repeated  indices  are  summed.  The  density  is 
subject  to  the  normalization  condition, 


yy*v(x1,x2)dx1dx2  = 1. 


1 2 

het  p(x  , x , t)  be  defined  by 


p (x1 , x2 , t ) = Prfx^t)  si  x1  (0)  = x1,  x2  (0)  = x2} 


Then  p will  satisfy  the  backward  equation  (Ludwig,  1975) 


pt  = + (bi  + °i,pi 


i 


subject  to  boundary  and  initial  conditions 


p (6 , x , t)  = 1 


lim  , 1 2 . . n 

1 P(x  , x , t)  = 0 

X -*«> 


ptx1,  X2,  0)  = 


0 x^>6 


1 otherwise 

Finally,  we  consider  the  expected  time  to  cross  a given  line 
T(6;x^,x2)  = E{t:x1(t)  £ SjX^ts)  > 6,s  < tlx^O)  = 

xQ,  x^(0)  = x^,  x (t)  eventually  crosses  6} 


(40) 


(41) 


Then  T satisfies 


/ 1 2, 
■P(x0'xo) 


| + (b1  + c1)Ti 


(42) 


1 2 

where  p(x  , x ) is  the  steady  state  solution  of  (38)  (e.g., 

o o 

Gihman  and  Skorohod (1972) ) . In  the  case  that  we  use  the  steady 


state  (e=0)  assumption 
as  will  be  seen. 


equations  (36),  (39),  (42)  simplify, 
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SECTION  3 


A 

A TWO-TIMING  METHOD  FOR  n,q 


We  now  return  to  the  deterministic  equations  (10), (12)  and 

analyze  them.  In  section  4 , 5 we  consider  the  stochastic  equations. 

Equations  (10) , (12)  represent  a singular  perturbation  problem.  It 

is  clear  that  for  times  x » e,  q has  achieved  its  steady  state  q 
A 

so  that  q = 0.  Our  goal  in  this  section  is  to  calculate  the  rate 
A 

at  which  q O.  We  introduce  an  "initial  layer"  time  scale 


S = - 

e 


Then  equations  (10) , (12)  become 


= e{n(l  - n)  - + anq}  n(0)  = Nc 


dn 

dS 


A 

da 

dS 


= - (n  + X)q 


A 


q (0)  = Q 


We  seek  a solution  of  the  form 


n(S)  = nQ(S)  + en1(S)  + e n2(S)  + 


q(S)  = qQ(S)  + eq^  (S)  + e + 


1 


After  the  series  (46,  47),  are  substituted  into  (44,  45),  terms 
are  collected  according  to  powers  of  e . We  use  the  expansion 


* • Jk  (•-%*-) 


The  0(1)  terms  indicate  that 


(48) 


dn 

o 

~dS 


= O 


no  = No 


(49) 


A 

dq 


o 

dS 


- (n  + X ) q 
o ^o 


-<No  + X)<3o 


Thus, 


qo(s) 


Qoexp[-(NQ  + A)S] 


(50) 


(51) 


To  leading  order  in  e,  q(S)  relaxes  towards  O with  a relaxation 
time 


N + A 
o 


(52) 


Corrections  to  the  above  relaxation  time  can  be  obtained  by 
considering  further  terms  in  the  series.  The  0(e)  terms  indicate 
that. 
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dn. 


dS 


aAN  A 

N -N  - - + aN  q 

o o A+N  o^o 


A 

dqj^ 

~dS 


= -q  (N  + A)  - q n 


* A 1 

O 1 


These  equations  do  not  have  a simple  solution. 

Next  we  consider  the  "outer"  equations,  in  the  time  variable  x 


We  seek  solutions  of  (10,  12)  of  the  form 


n (x ) = £ e^n.(T) 


j=0 


q(T)  = £ e-’q-tt) 


j = 0 


The  0(1)  terms  indicate  that 


dn 


dx 


- -2  “>no 

n - n 

o o 


+ an  q 


A+n 


o o 


0 = 


-q  (n  + A) 
^o  o 


(53) 


Hence,  we  obtain  the  steady  state  result  of  Clark  and  Mangel  (1977) 


dno  - -2 

<rr~  = no  " no 


aAn 


A+n 
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(54) 


(55) 


(56) 


(57) 


(58) 


(59) 


l 


No  simple  solution  for  nQ(T)  is  available  in  terms  of  elementary 
functions.  Consequently,  the  usual  "matching"  procedures  of  singular 
perturbation  theory  (e.g.,  Lin  and  Segel  (1977))  are  not  very  helpful 
(Simple  considerations  show  that  the  matching  procedure  must  yield 

vo)  - v> 

We  now  consider  equation  (59)  for  the  special  case  that  a > 1 
and  A > XQ,  the  bifurcation  value  of  X.  In  this  case,  the 
origin  is  the  only  steady  state.  Consequently,  the  population 
is  driven  to  extinction.  Let 


T(6,  Nq)  = min{t :nQ (t) ^6 , nQ(S)>6  S<t,  nQ(0)  = Nq} 


The  interpretation  of  T(6,  Nq)  is  as  an  "extinction  time",  if 
6 is  small.  We  will  now  construct  a second  perturbation  expansion, 
for  large  X,  in  order  to  calculate  T(6,  Nq) . 

Equation  (59)  can  be  inverted  to  give: 


dT 

dn_ 


X + n 


n (1  - n ) (A  + n ) - a An 
o o o o 


T(Nq)  = 0 


We  introduce  f = 1/A  as  a second  small  parameter,  so  that 


dT 

d*o 


1 + rn 

o 

n (1  - n ) (1  + Tn  ) - an 
o o o o 


(60) 


(61) 


(62) 
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We  seek  a solution  of  (62)  in  the  form 


T = £ r3Tj(no) 


(63) 


The  leading  terms  satisfy 


dT 


dn  n (1  - n ) - cm 

o o o o 


T (N)  = 0 
o o 


(64) 


dT 


1 _ 


n 


n2*1  - V 
o 


_ - _ - _ — Y 

dn  n ( 1 - n ) - an  [n  (1  “ *0  ~ an;J 

o o o o o o o 


T1(No)=0  (65) 


These  equations  can  be  easily  solved.  Thus,  we  find: 

T(s'  *V  - A [ln(jr)  - ln(i^s;)] 

+ “X  | (0-1)  + 6 " (0-1)  + Nq  | + 0 (^2 ) 


(66) 


In  summary,  the  result  (66)  was  obtained  by  using  two  perturbation 
procedures.  First,  we  used  a singular  perturbation  procedure  to 
calculate  the  relaxation  time  tr  and  to  obtain  (59).  Second,  a 
regular  perturbation  procedure  was  used  to  calculate  T(6,  Nq) . 
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To  check  the  accuracy  of  the  steady  state  assumption,  T(6,  Nq) 
calculated  from  (66)  was  compared  with  the  exact  solution  of  (10,13), 
by  a Runge-Kutta  method.  We  chose  q(0)  =1,  e = 1,  6 = .1.  In 
table  1,  some  results  of  the  calculation  are  presented.  Since 
e = 1,  in  principle, the  use  of  the  steady  state  assumption  is  not 
warranted.  However,  the  results  in  table  1 indicate  that  the 
steady  state  assumption  provides  a relatively  good  description  of 
the  system,  even  for  e order  1. 


n 


SECTION  4 

EXTINCTION  IN  A RANDOM  ENVIRONMENT:  STEADY  STATE  CASE 

We  now  consider  some  of  the  stochastic  problems  presented 

in  section  2.3,  assuming  that  the  steady  state  assumption  (e  = O) 
A 

about  q holds.  Hence,  the  stochastic  equations  (31,  32)  are 


replaced  by. 


dn 


= £n(l  - n)  - dt+/aa  (n) 


dW 


where  a = 1/N  and 

, x 2 

a (n)  = n + n J 


gin 

A+n 


When  a = 0 (i.e.,  1/N  = 0)  we  obtain  the  deterministic  kinetic 
equation  (59)  from  (67).  This  is  in  accord  with  intuition:  as 
the  sample  size  becomes  larger,  fluctuation  effects  should  dis- 
appear. 

The  flow  of  the  deterministic  equation  is  sketched  below 
(these  pictures  are  transcriptions  of  the  results  in  section  2) : 


n 


nT 


4 


n 


a < 1,  all  A 

u > 1 , A < A 


H 


a > 1 , A > A 


B 


B 


(67) 


(68) 


In  order  to  simplify  calculations,  we  shall  approximate  a (n)  by  n. 

Since  the  stochastic  effects  are  most  interesting  when  n is  small, 
this  approximation  (which  greatly  simplifies  calculations)  seems  reason- 
able. 
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THE  DENSITY  v(n,  t) 


The  forward  equation  for  v(n,  t)  (35)  becomes, 


v. 


% (llVlnn 

2 nn 


c- 


(1-n)  - 


aXn 

A+n 


a, 

4*v 

-I  n 


The  first  question  we  wish  to  consider  is  the  calculation  of  the 
stationary  distribution,  v(n)  . If  vfc  - 0 , (68)  becomes  an 

ordinary  differential  equation  that  can  be  integrated  to  yield. 


v (n) 


In  n 


The  constant  c is  determined  by  normalization: 


c 


-1 


L J 

The  probability  of  eventually  finding  the  population  at  a level 
less  than  6 is. 


P 


6 


v (n)dn. 


Natural  questions  concern  the  location  of  concentration  points  of 
v(n)  and  the  behavior  of  v(n)  for  A large.  The  density  v(n) 


(68) 

(69) 

(69a) 

(70) 

(71) 
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will  be  concentrated  near  the  minima  of  the  function  <}>  (n)  . It 
is  easy  to  show  that  these  minima  correspond  to  the  steady  states 
of  the  deterministic  kinetic  equation. 

When  X is  large,  we  expand  ln(l  + n/X)  in  its  Taylor 


series  and  obtain. 


v(n)  = c 


exp  r |n  - £ - «»(?  - £r) + - a4 »j 


Thus,  for  small  a and  large  X: 


v(n)  ~ c 


expf o~  in(1  _ a)  " ^ ' r)}] 


In  this  case,  the  extremum  of  <(>  occurs  at 


n*  = 


The  result  (76)  indicates  the  following.  If  a < 1,  then  n*  > 0 
always,  so  that  the  density  is  concentrated  near  n*.  When  a > 1, 
n*  > 0 if  X is  small  enough.  In  particular,  n*  > 0 if 


1 - Y < 0 

A 


i . e . , X < a 


Equation  (76)  can  be  viewed  as  a stochastic  bifurcation  result.  We 
note  that,  for  a > 1,  n*  corresponds  to  a minimum  of  v(n). 
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MEAN  TIME  TO  CROSS  n = 6 

We  now  turn  to  the  calculation  of  the  mean  time  to  cross  the 
level  n = 6,  conditioned  on  n(0)  = Nq  . This  function  T(6,  Nq) 
satisfies 

Xfnn  + - X5T+  ?]fn  = -?<"> 

In  (77),  p(n)  = Pr  {n(t)  eventually  crosses  6|n(0)  = n} . 

The  function  T satisfies  the  boundary  conditions: 


T (6 , 6)  = 0 


lim  Tn (6  ,n)  = 0 

n->oo 


(78) 


The  function  p(n)  satisfies: 


°D_  p + 
2 ^nn 


a \ ctXn  , a _ _ „ 

n(1-n>  ' X+n  * 1 |pn  - 0 


(79) 


subject  to 

P (6)  = 1 


lim  & p (n)  ->0 
n—m 


(80) 


Both  the  functions  p(n)  and  T(6,  n)  can  be  calculated  explicitly 
by  two  integrations: 

2 


p (n) 


/ exp  -|jn'  - ^ 


- aXln (1+ 


| ln  n'[ 


dn’ 


(81) 


j f exp[-ljn’  - ^ dnl 
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1 


In  many  cases  of  interest,  p(n)^  1 for  almost  all  values  of  n. 
The  solution  of  the  problem  posed  by  (77,78)  is, 


^ r ( 2 

(n)  = £ J exp  ~^js  “ ^2  _otA  - -j  0n  s| 


(82) 


°°  2 

jf  exp[f  jy  - \ _aA  **  (1+x>  ■ xf  ^ yj 


p (y)dydS 


The  behavior  of  p(n),  T(6,  n)  is  sketched  in  figure  4,  for  a 
number  of  cases.  Equations  (81,  82)  allow  the  explicit  calculation 
of  the  mean  time  to  extinction  in  the  stochastic  case. 
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SECTION  5 


EXTINCTION  IN  A RANDOM  ENVIRONMENT, 

WITHOUT  THE  STEADY  STATE  ASSUMPTION 

A 

In  two  dimensions  the  [ (n,q)  or  (n,q)  phase  space]  no 
exact  solutions  of  the  stochastic  problems  are  possible.  If  the 
intensity  of  the  noise  is  small  (i.e.,  1/N,  1/Q  large),  then 
approximate  theories  are  available  (Ludwig  (1975),  Mangel  (1977)). 

In  this  section,  we  will  sketch  how  those  theories  apply  to  the 
aggregation  problem.  We  wish  to  calculate:  1)  the  probability 
that  n(t)  ever  crosses  the  line  n = 6,  conditioned  on  initial 
value  of  n,  and  2)  the  mean  time  to  cross  the  line  n = 6.  We 
differentiate  two  cases,  according  to  the  values  of  a. 

a < 1;  THE  RAY  METHOD  (LUDWIG , 1975) 

When  a < 1 (figure  1) , the  non-trivial  steady  state  is 
always  stable.  Hence,  we  are  interested  in  the  exit  from  the  domain 
n > 5.  The  density  of  (n,q) , v(x,n,q),  satisfies  the  forward 
equation 

vT  = | (ai^v)ij  - ( (b1  + aci)v)i  (83) 


(84) 


with  a(n,q),  c(n,q)  given  by  (33),  (36)  respectively  and 

b^  = , b2  = in  (10,12).  We  assume  that  e = 1. 

dx  dx 

Ludwig  (1975)  has  shown  that  the  ray  method  (Cohen  and  Lewis,  1967) 
applies  here.  A solution  of  (83)  of  the  form, 

v(x,n,q)  = expj^i^MiJ  £ Z.oj 
is  sought. 


L 


I 

I 
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After  derivatives  are  evaluated,  terms  are  collected  according  to 
powers  of  a.  The  leading  term  vanishes  if 


<b  . d>  • 


o 


(85) 


To  first  order,  Zq  satisfies  an  ordinary  differential 
equation  on  the  characteristics  of  (85) . 

Equation  (85)  can  be  solved  by  the  method  of  characteristics,  once 
initial  data  are  given.  Ludwig  has  shown  how  this  may  be  accomplished. 
The  characteristics  of  the  system  (85)  are  called  rays. 

Suppose  that  v ~ 0 on  the  boundary  and  that  we  assume  for  v. 


v (t  , n , q) 


2-r  e v(n,q) 

k=0 


(86) 


Then  the  expected  time  to  hit  the  boundary  is  asymptotic  to  1/XQ 
(Ludwig,  1975)  . We  are  led  to  the  eigenvalue  problem 


-X  = % (aljv) . . - ( (b1  + oc1)v).  (87) 

o 2 lj  l 

Ludwig  has  shown  how  Miller's  method  (Miller,  1962)  may  be  generalized 
to  calculate  XQ.  Thus,  in  principle  the  problem  is  completely 
solved.  Further  details  and  calculations  are  in  (Ludwig,  1975) . 

a a 1:  THE  GENERALIZED  RAY  METHOD  (MANGEL,  1977) 

When  a a 1,  there  may  be  two  non-zero  steady  states  (X  < XR) , 
degenerate  steady  state  (X  = Xg)  , or  only  the  trivial 


one 


steady  state  (X  > Xn)  (figure  2) . The  ray  method  will  break  down 
in  this  case,  but  a generalization  of  the  ray  method  does  work 
(Mangel,  1977).  Here  we  sketch  the  technique.  The  details  can  be 
found  elsewhere.  We  consider  T(6;  n,q) , which  satisfies 


-1 


+ b1T.  + 
i 


c 


(88) 


Asymptotic  analysis  of  a one  dimensional  problem  (for  a small) 
indicates  that  a formal  solution  of  (88)  is 


T = g (x)  B (i^/a3^3  , a /a2//3  , l/cs1//3,  n) 


+ e2/3h(x)B'  (ij>/a1/3,  a/a2/3,  l/a1/3,  n) 


(89) 


+ ek(x) 


In  (89),  B (Z , a, f ^ ) satisfies  the  differential  equation 


d2B 


dZ 


,z2  - “’at  - q * qz  zo  5 z 


(90) 


The  functions,  g(x),  h(x),  and  k(x)  should  also  be  expanded  in 
asymptotic  series  of  the  form. 


g(x)  = ^ engn(x) 


k(x)  = S enkn(x) 


h (x ) = £enhn(x)  (91) 
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where  the  superscript  n is  an  index  on  gn,hn,  and  kn.  After 
derivatives  are  evaluated,  and  substituted  into  (88)  (using  equation 
(90)  to  eliminate  B"  and  B ' ' * ) terms  are  collected  according  to 
powers  of  e.  The  leading  coefficients  vanish  if 

0(y"2/3B):  - -2-  ('f'2  - ct)  = 0 (92) 

O ( Y°B) : b1gi  = 0 (93) 

i j 

O ( Y°)  : b1ki  + ^2“  g^-j  (-1  + .pn)  = "I  (94) 


Equation  (92)  is  a generalized  Hamilton-Jacobi  equation.  It  can  be 

solved  by  the  method  of  characteristics.  The  parameter  a is 

2 

determined  by  requiring  i(i  = a at  the  points  where  b(x)  vanishes. 

2 

Namely,  at  the  steady  states  Pq,  P^,  iJj  = a.  The  value  of  a 
must  be  determined  by  a numerical  iterative  procedure  (Mangel,  1977  ). 
At  \ = A0,  so  that  Pn  and  P,  coalesce  , the  value  of  a = 0. 
Equation  (93)  indicates  that  g is  a constant.  It's  value  is 
determined  as  follows.  At  PQ,  P^,  equation  (94)  becomes 


l j 

g— 5—  iK<K(-l  + i^n)  | = -l 

1 J D 

*k 

For  k = 0 , 1 , we  obtain  two  equations  for  g and  n • When 
\ = A , it  is  possible  to  show  that  n = 0 and  (95)  provides 
one  equation  for  g.  (Mangel,  1977)  . 


(95) 
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— ■■ 


We  give  data  k = 0 on  n = 6.  Then  (94)  can  be  solved  by  the 

method  of  characteristics.  We  let  4»  = on  n = 6,  set  = ZQ 

in  (90)  and  B(Z  ) = B ' ( ZQ ) = 0.  Then  T = 0 on  n = 6.  This 

2/3 

construction  will  yield  values  for  T which  are  0(o  ) different 

from  the  true  values.  Further  details  can  be  found  in  Mangel  (1977  ) . 
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APPENDIX  A 


THE  GENERAL  AGGREGATION  MODEL 

In  this  appendix,  we  prove  a theorem  that  applies  to  the 
general  set  of  aggregation  equations 

N = gN  - af(N  ,Q) 

Q = af(N,Q)  - XQ 

with  N,  Q defined  as  in  the  body  of  the  paper. 

Let  Q*  {N,  « , X ) be  the  solution  of 

uf(N,Q)  = XQ 

such  that  the  following  hypotheses  are  satisfied: 


a. ) 

lim 

Q*  (N, 

a , X) 

V 

c 

o 

ii 

00 

N-»°° 

b.) 

lim 

N-.00 

of  (N, 

Q*  (N, 

e 

ii 

o 

i—* 

o' 

>• 

A 

8 

for 

all  a,  X 

c. ) 

lim 

X-.°° 

Of  (N, 

Q*  (N, 

a,  X))  - 

C2  (a , N)  < 00 

for 

all  a,  N 

Then  there  exist  values  a*,  X*  such  that  the  tuna  steady  state 
level  exhibits  a bifurcation.  Namely,  for  a > a*  and  X < X* 


A-l 


(A-l) 

(A-2) 

(A-3 ) 


has  two  solutions:  N = 0 and  N = N*  , where  n*  is  determined  by 

j|y<N,  AMI  -0  • (A-9) 

N* 

Then  for  fixed  a,  a > a*,  as  X increases  through  X*  , we 
observe  a bifurcation  (figure  5) . 


APPENDIX  B 


A KINETIC  MODEL  FOR  FLUCTUATIONS 

Clark  and  Mangel  considered  a kinetic  model  that  yields  the 
same  behavior  as  equations  (11,  12) . This  model  provides  a natural 
way  to  treat  fluctuations.  We  let  T(t)  be  the  number  of  "core" 
schools  in  the  ocean  and  consider  the  kinetic  scheme. 


bEX 

C^  -•  ° harvest  + K 
bEX 

C2  -•  0 harvest  + K 

where  C-^,  C2  are,  respectively,  surface  complexes  with  1 and 
2 core  schools.  The  kinetic  equations  for  the  scheme  (B-l) , derived 
by  a "law  of  mass  action"  approach,  are: 

gT  - alTK  - a2TC;L  + B^  + B^ 

gR  - n1TK  + 0^  + bEX  (Cj  + C2) 


dT 

dt 


dK 

dt 


(B-l) 


(B-2 ) 


(B-3) 


i 

i 


1 


(B-4) 


= o^KT  - &1C1  + e2C2  - a2C1T  - bEX^ 

= a2ClT  - B2C2  - bEXoCl  (B-5) 

The  analysis  of  (B.2-B.5)  yields  results  analogous  to  (10-12); 
if  we  set  B^  = 62  = 0 (Clark  and  Mangel,  Appendix  B) . 

For  the  purpose  of  stochastic  modeling,  however,  the  scheme 
(B-l)  has  a number  of  advantages.  In  particular,  no  re-interpretation 
of  the  kinetic  equation  is  needed,  since  the  elementary  events  (B-l) 
are  clearly  a birth  and  death  process.  Consider  a time  increment 
t -*■  t + At  , with  the  values  of  T(t),  C^(t),  C2(t)  and  K ( t ) 
known.  Then  in  the  increment  At,  we  assume  that  the  following 
transition  probabilities  exist: 


Pr{  AT  = 1,  AK  = 1,  A^  = -1,  AC2  = 0}  = B-jCjAt  + o(At)  (B-6) 

Pr{  AT  = -1,  AK  = -1,  A^  = 1,  AC2  = 0}  = c^TKAt  + o(At)  (B-7) 

Pr{ AT  = -1,  AK  = 0,  A^  = -1,  AC2  = 1}  = a2TC1At  + o (At)  (B-8) 

Pr{ AT  =1,  AK  = 0,  LC1  = 1,  AC2  = -1}  = B^At  + o(At)  (B-9) 

Pr{ At  =0,  Ak  = 0,  A^  = -1,  AC2  = 0}  = bEX^At  + o(At)  (B-10) 


B-2 


( B — 1 1 


Pr{ AT  = 0,  AK  = 0,  ACj.  = 0,  AC2  = -1}  = bEX^At  + o(At) 

Pr{all  other  transitions}  = o(At), 

where  AT  = T(t  + At)  - T(t),  etc. 

If  we  use  a quasi-steady  state  approach,  we  replace  C^,  C2>  and  K 
by  the  steady  state  values  (Clark  and  Mangel (1977 )) : 


K - % 

(B-13 ) 

c _“2C1T 

( B— 1 4 ) 

^2  S-  + bEX 

2 o 

a 

jTK 

(B-15) 

^1  "" 

6i  + a2 

T + bEX  a262T  \ 

° l»2*  bEXo  ! 

Then  the  transition 

probabilities  are  functions  of  T(t)  only. 

We  obtain 

Pr{ AT  = 1}  = 

At{gT  + B1C1  + e2C2>  + o (At) 

(B-16) 

= 

X (T)  At  + o (At) 

(B-17 ) 

Pr{ AT  = -1}  = 

AttcxjKT  + OjCjT} 

(B-18 ) 

Vi  (T)  At  + o ( At) 


with  , C2>  and  K given  by  (B. 13-15).  It  is  easy  to  show  that 


Ito  re  E,4Tl 


= X (T)  - V (T) 


is  the  same  as  the  law  of  mass  action  result.  We  can,  however, 
immediately  calculate  the  incremental  variance  needed  for  the 
diffusion  approximation: 


a(T>  = it^o  a^e{(at)2} 


= X (T)  + M (T) 


(B-19) 


(B-20) 


(B-21) 
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14.4 
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14.1 

4.21 

4.91 

14.2 

(c)  D < 0 

The  (n,q)  phase  plane  for  the  case  a>l.  In 
this  case,  the  deterministic  system  is  of  the  marginal  type 
(Mangel,  1977).  a)  For  D>0  (X<XB)  PQ  is  unstable.  (b  When 
D = 0 (X  = XB ) PQ  and  coalesce  to  form  a state  of 
marginal  stability.  c)  For  D<0  (X>XB)  PQ  and  Pi  have 
annihilated  each  other.  The^only  steady  state  is  the 
origin;  it  is  stable.  The  n = 0 is  really  more 
parabolic  than  shown. 
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